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Abstract 

We investigate the origin of the deviations from the classical Darcy law by 
numerical simulation of the Navier-Stokes equations in two-dimensional dis- 
ordered porous media. We apply the Forchheimer equation as a phenomeno- 
logical model to correlate the variations of the friction factor for different 
porosities and flow conditions. At sufficiently high Reynolds numbers, when 
inertia becomes relevant, we observe a transition from linear to non-linear be- 
havior which is typical of experiments. We find that such a transition can be 
understood and statistically characterized in terms of the spatial distribution 

of kinetic energy in the system. 
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A standard approach in the investigation of single-phase fluid flow in microscopically 
disordered and macroscopically homogeneous porous media is to characterize the system in 
terms of Darcy's law which assumes that a global index, the permeability k, relates 

the average fluid velocity V through the pores with the pressure drop AP measured across 
the system, 

Here L is the length of the sample in the flow direction and /i is the viscosity of the fluid. 
However, in order to understand the interplay between porous structure and fluid flow, 
it is necessary to examine local aspects of the pore space morphology and relate them to 
the relevant mechanisms of momentum transfer (viscous and inertial forces). This has been 
accomplished in previous studies [^-[T0[| where computational simulations based on a detailed 



description of the pore space have been quite successful in predicting permeability coefficients 
and validating well-known relations on real porous materials. 

In spite of its great applicability, the concept of permeability as a global index for flow, 
should be restricted to viscous flow conditions or, more precisely, to small values of the 
Reynolds number. Unlike the sudden transition from laminar to turbulent flow in pipes 
and channels where there is a critical Reynolds number value separating the two regimes, 
experimental studies on flow through porous media have shown that the passage from linear 
(Darcy's law) to nonlinear behavior is more likely to be gradual (see Dullien and references 
therein). It has then been argued [Q and conffimed by numerical simulations that 



the contribution of inertia to the flow in the pore space should also be examined in the 
framework of the laminar flow regime before assuming that fully developed turbulence effects 
are present and relevant to momentum transport. Here we show by direct simulation of the 
Navier-Stokes equations that the departure from Darcy's law in flow through high porosity 
percolation structures (e > ec, when ec is the critical percolation porosity) and at sufficiently 
high Reynolds numbers can also be explained in terms of the inertial contribution to the 
laminar fluid flow through the void space. The calculations we perform do not apply for 



unstable or turbulent Reynolds conditions. We then demonstrate that it is possible to 
statistically characterize the transition from linear to nonlinear behavior in terms of the 
distribution of kinetic energy. This allows us to elucidate certain features of the fluid flow 
phenomenon in irregular geometries that have not been studied before. 

Our model for the pore connectivity is based on the general picture of site percolation 
disorder. Square obstacles are randomly removed from a 64x64 square lattice until a porous 
space with a prescribed void fraction e is generated. The mathematical description for the 
detailed fluid mechanics in the interstitial pore space is based on the assumptions that we 
have steady state flow in isothermal conditions and the fluid is continuum, Newtonian and 
incompressible. Thus, the continuity and Navier-Stokes equations reduce to 

V-u = , (2) 

p u ■ Vu = -Vp + fi V^u , (3) 

where p is the density of the fluid and u and p are the local velocity and pressure flelds, 
respectively. We use the nonslip boundary condition at the whole of the solid-fluid interface. 
End effects of the flow fleld established inside the pore structure (particularly signiflcant at 
high Reynolds conditions) are minimized by attaching an inlet and an outlet to two oppo- 
site faces. At the inlet a constant inflow velocity in the normal direction to the boundary is 
specified, whereas at the outlet the rate of velocity change is assumed to be zero (gradient- 
less boundary condition). Instead of periodic boundary conditions, we close the remaining 
two faces of the system with two additional columns of obstacles. This insulating condition 
reproduces more closely the experimental setup usually adopted with real rocks and per- 
meameters. The Reynolds number is defined here as Re = pVdp/fi where dp is the grain 
diameter . For a given realization of the pore geometry and a fixed Re, the local velocity 
and pressure fields in the fiuid phase are numerically obtained through discretization (see |T^ 



for numerical details) by means of the control volume finite-difference technique [p!^ , p!5[] . Fi- 
nally, from the area-averaged pressures at the inlet and outlet positions, the overall pressure 
drop can be readily calculated. 



The classical approach to macroscopically characterize the effect of inertia on flow 
through real porous media is to use the Forchheimer equation [|1],0, 

A p 

- = a^V + (3pV'^. (4) 



The coefficient a corresponds to the reciprocal permeability of the porous material and 13 is 
usually called the "inertial parameter" . Both a and jd should depend on the porosity e of the 
porous material. At sufficiently low velocities, Eq. (4) reduces to Darcy's law, Eq. (1). The 
term jSpV"^ can be interpreted as a second order correction to account for the contribution of 
inertial forces in fluid flow. Equation (4) is not a purely empirical expression, since it can be 
derived by an appropriate average of the Navier-Stokes equation for one-dimensional, steady 
incompressible laminar flow of a Newtonian fluid in a rigid porous medium Rearranging 
(I) in the form, 

where / = —AP/L/SpV"^ and Re' = jSpV/ap, we obtain a friction factor- Reynolds number 
type of correlation which is presumably "universal." Equation (5) has been extensively and 
successfully used to correlate experimental data from a large variety of porous materials 
and a broad range of flow conditions . Certainly, a better representation for experimental 
data in the non-Darcy flow regime can be obtained with the addition to the Forchheimer 



equation of third order corrections in the velocity P Jll|Jl^ . The theoretical basis for this 
type of correction term, however, is still controversial. 

Figure 1 shows the results of our flow simulations in terms of the Forchheimer variables / 
and Re' for three different values of lattice porosity (e = 0.7, 0.8 and 0.9). After computing 
and averaging the overall pressure drops for all realizations at different values of e and Re, 
we fit the results with Eq. (4) to estimate the coefficients a and (3 and calculate / and Re'. 
In agreement with real flow experiments, we observe a transition from linear (Darcy's law) 
to nonlinear flow. Moreover, the point of departure from linear to nonlinear behavior in the 
range 10~^ < Re' < 10~^ is consistent with previous experimental observations. However, 



in the region 6.2 x 10~^ < Re' < 1.8, the Forchheimer equation generally overestimates the 
computed values of the friction fraction. In addition, for a fixed value of Re', the variability 
in this transition region is sufficient to suggest a dependence of the type / = /(-Re', e) |T6| . 

The flow distribution in two-dimensional incompressible systems can be conveniently 
described in terms of the stream function ip [|T7|. Figure 2a shows the contour plot of ip 



for a typical realization of a highly porous void space (e = 0.9) subjected to low Reynolds 
conditions. Re = 0.0156. In spite of the well-connected pathways available for flow at this 
large porosity value, the predominant viscous forces in the momentum transport through 
the complex void geometry generates well defined "preferential channels" of fiuid fiow. As 
shown in Fig. 2b, the situation is quite different at high Re, where the degree of channeling 
is clearly less intense than in Fig. 2a. In the case of Fig. 2b, due to the relevant contribution 
of inertial forces to the fiow, the distribution of streamlines along the direction orthogonal 
to the main fiux y becomes more homogeneous. 

The channeling effect can be statistically quantified in terms of the spatial distribution 
of kinetic energy in the flowing system. In analogy with previous work on localization of 



vibrational modes in harmonic chains |T3], we deflne a "participation" number vr 



where n is the total number of fluid cells in the numerical grid enclosing the physical pore 
space, Qi = ei/J2j=i ^i; ^ (""i^ + "^j^) kinetic energy associated with each individual 

fluid cell, and Ui and Vi are the components of the velocity vector at cell i in the x and 
y directions, respectively. From the definition Eq. (6), vr = 1 indicates a limiting state of 
equal partition of kinetic energy (gj = 1/n, Vi). On the other hand, a sufficiently large 
system {n oo) exhibiting strong channeling effects should correspond to a "localized" 
fiow field, TT ^ [|19|. We calculate the function tt for 10 pore space realizations generated 
with e = 0.9 at different Re. Due to convergence difficulties and computational limitations 
on the resolution of the numerical grid, we restrict our calculations to Re < 15.6. As shown 
in Fig. 3, the participation number remains constant, vr ^ 0.37, for low Re up to a transition 



point at about Re ^ 0.3. Above this point, the flow becomes gradually less localized (tt 
increases) as Re increases. This transition reflects the onset of inertial effects in the flow, 
and the significant changes in tt above the transition point indicate the sensitivity of the 
system to these nonlinearities. The large error bars at low Re indicate that tt is sensitive 
to structural disorder if the viscous forces are effectively generating preferential channels in 
the fiow. 

The difference between our results at low and high Re can be better understood if we 
remember that viscous effects extend a long way at low Reynolds conditions, so that distant 
boundaries may have a large effect on the streamlines. It is then interesting to visualize 
the distortions in the local velocity field when inertial forces become important compared to 
viscous forces. In Fig. 4 we show the profiles of the velocity magnitude at different positions 
along the main fiow direction a; in a typical realization of the porous media. At low Re 
(Fig. 4a), the fiuctuations in the velocity field are essentially smooth in shape, with peaks 
that closely correspond to the variations in the local porosity. At high Re (Fig. 4b), the 
situation becomes quite different. Due to inertia, the effect on the fiow field of the disorder 
in the local pore geometry tends to propagate further the fiuctuations in the x direction. We 
can follow in Fig. 4b the changes in shape of the velocity magnitude at different x positions. 
If there is an available straight void space pathway for fiuid fiow, a peak generated at a smaller 
X can persist further up in the next profiles located at larger x values. As a consequence, 
the velocity profiles at large Re are more rough than those at low Re. 

In summary, to characterize the infiuence of inertial forces on the flow of a single fiuid in 
porous structures, we demonstrate that incipient deviations from Darcy's law observed in 
several experiments can be modeled in the laminar regime of fluid flow, without including 
turbulence effects. The results of our simulations agree with numerous experimental data 
which display a gradual transition at high Re from linear to nonlinear flow in the pore space. 
Moreover, we show that this flow transition can be characterized in terms of the partition 
of kinetic energy in the fiuid phase. Namely, the fiow at low Re is more "localized" due 
to channeling effects than the flow at high Re conditions. Finally, our calculations with 
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the Navier-Stokes equations indicate that the Forchheimer model should be valid for low 
Re and also for a limited range of high Re numbers, even when inertial nonlinearities can 
significantly affect the momentum transport at the pore scale. However, the magnitude of 
the deviations we find at the transition from Darcy to non-Darcy flow suggests a nonuniversal 
behavior of the friction factor / with the porosity e within this particular region. 

We thank CNPq, FUNCAP and NSF for support. We also thank three referees for their 
constructive criticisms. 
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FIGURES 

FIG. 1. Logarithmic plot showing the dependence of the generahzed friction factor / on the 
modified Reynolds number Re' . The soUd hnc is the best fit to tlic Forchhcimcr equation (4), while 
the dashed line is the best fit to Darcy's law of the data at low Re' . Please note that three pairs 
of a and (3 parameters have been estimated for each of the three porosity values (e = 0.7, 0.8 and 
0.9). These simulations have been performed with up to five lattice realizations. The error bars 
are smaller than the symbols. 

FIG. 2. (a) Contour plot of the stream function ij) for low Reynolds number conditions 
{Re = 0.0156). (b) Same as in (a), but for a Reynolds number 1000 times larger [Re = 15.6). In 
both plots, the values of ip change by a constant increment between consecutive streamlines. 

FIG. 3. Dependence of the participation number tt on the Reynolds number Re (e = 0.9). 
These simulations have been performed with ten lattice realizations. 

FIG. 4. (a) Profiles of the velocity magnitude at different positions in direction x for low 
Reynolds number conditions {Re = 0.0156). The realization is the same as in Fig. 2, but the 
obstacles have been removed for better visualization, (b) Same as in (a), but for a Reynolds number 
1000 times larger {Re = 15.6). The velocity magnitudes in (a) and (b) have been normalized by 
the maximum value calculated for each system. 
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